1 Data import

The data is included within this reproducible report and can be downloaded here (object 'rbt_public_1') and here (object 'rbt_public_0').

1.1 Data wrangling

# wrangling
rbt_public_l <- pivot_longer(rbt_public_1, 
                             cols = (abs1_tsm_1:abs2_tch_5), 
                             names_to = "variable", 
                             values_to = "value") %>%
    mutate(treat = case_when(                    # create treatment variable
        str_detect(variable, "abs1_") ~ toupper(treat1),
        str_detect(variable, "abs2_") ~ toupper(treat2)),
        variable = str_sub(variable, 6, -1)) # delete title page substring

# pivot them into wide format
rbt_public_w <- pivot_wider(rbt_public_l,
                            names_from = variable,
                            values_from = value,
                            id_cols = c(session, treat))

# invert trust items & build scales
inv7fct <- function(x) (8-as.numeric(x))
data_rbt <- rbt_public_w %>%
    mutate_at(vars(tru_exp_1:tru_ben_4),
              list(~inv7fct(.))) %>% # recoding 1-->7, 2-->6, ...
    rename(exp_1 = tru_exp_1, 
           exp_2 = tru_exp_2, 
           exp_3 = tru_exp_3, 
           exp_4 = tru_exp_4, 
           exp_5 = tru_exp_5, 
           exp_6 = tru_exp_6,
           int_1 = tru_int_1, 
           int_2 = tru_int_2, 
           int_3 = tru_int_3, 
           int_4 = tru_int_4,
           ben_1 = tru_ben_1, 
           ben_2 = tru_ben_2, 
           ben_3 = tru_ben_3, 
           ben_4 = tru_ben_4) %>% 
    group_by(session) %>% 
    mutate(meas_rep = 1:n()) %>% 
    ungroup() %>% 
    full_join(., rbt_public_0)

## Joining, by = "session"

# data frames with sum scales
data_scales <- data_rbt%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
         TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))

data_scales_lables <- data_rbt%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Treatment = fct_recode(Treatment,
                                `Colored badges` = "CB",
                                `Control Condition` = "CC",
                                `Greyed out badges` = "GB"),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
        `Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), 
                                               na.rm = T))

data_scales_wide <- data_scales%>%
  select(Experitse, Integrity, Benevolence,
         TSM, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

data_scales_lables_wide <- data_scales_lables%>%
  select(Experitse, Integrity, Benevolence,
         `Topic specific multiplism`, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

# Data of treatment check
data_tch <- data_rbt %>% 
  select(starts_with("tch"), meas_rep, treat)

data_tch_n <- data_tch%>%
  mutate_all(function(x) ifelse(x == -999, NA, x)) %>% 
  filter(rowSums(is.na(data.frame(tch_1, tch_2, tch_3, tch_4))) != 4)

PID_list <- data_rbt$session %>% unique()

2 Sample

2.1 Sample size

length(PID_list)

## [1] 257

data_rbt %>%
  filter(meas_rep == 2) %>%
  select(session) %>% 
  distinct() %>% 
  nrow()

## [1] 257

2.2 Skipped after first measurement

sum(is.na(data_rbt%>%
            filter(meas_rep == 2)%>%
            pull(exp_1)))

## [1] 0

2.3 Demographics

data_rbt %>% 
  filter(meas_rep == 2)%>%
  mutate(education = as.factor(education),
         sex = as.factor(sex),
         age = as.factor(age)) %>% 
  select(age, education, sex) %>% 
  skim(.)

Data summary
Name	Piped data
Number of rows	257
Number of columns	3
_______________________
Column type frequency:
factor	3
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
age	1	FALSE	3	50: 108, 16: 78, 35: 71
education	1	FALSE	13	app: 61, L34: 56, L12: 26, L12: 19
sex	1	FALSE	2	f: 134, m: 123

2.3.1 Count on sex

data_rbt %>%
  filter(meas_rep == 2)%>%
  pull(sex) %>% 
  table(.)

## .
##   f   m 
## 134 123

2.3.2 Count on education

data_rbt %>%
  filter(meas_rep == 2)%>%
  pull(education) %>% 
  table(.)

## .
## app_11 app_12 app_13  app_5  L12_1  L12_2  L12_3  L12_4 L34_10  L34_6  L34_7 
##     11      3     61      6     19      9     26     18     14      8     18 
##  L34_8  L34_9 
##     56      8

2.3.3 Count on age

data_rbt %>%
  filter(meas_rep == 2)%>%
  pull(age) %>% 
  table(.)

## .
##  16  35  50 
##  78  71 108

3 Psychometric properties of the measurements

3.1 Descriptive overview

skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))

my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))

data_scales%>%
  my_skim(.)

Data summary
Name	Piped data
Number of rows	514
Number of columns	37
_______________________
Column type frequency:
character	4
factor	1
numeric	32
________________________
Group variables	None

Variable type: character

skim_variable	complete_rate	min	max	n_unique
session	1	64	64	257
treat	1	2	2	3
sex	1	1	1	2
education	1	5	6	13

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
Treatment	0	1	FALSE	3	GB: 183, CC: 170, CB: 161

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist	skewn	kurto	maxabszvalue
tsm_1	1	2.52	0.79	1	2.00	3.00	3.00	4	▂▇▁▇▂	-0.11	-0.43	1.92
tsm_2	1	2.72	0.82	1	2.00	3.00	3.00	4	▁▇▁▇▃	0.05	-0.74	2.11
tsm_3	1	2.56	0.83	1	2.00	3.00	3.00	4	▂▇▁▇▂	0.02	-0.58	1.88
tsm_4	1	2.82	0.82	1	2.00	3.00	3.00	4	▁▅▁▇▃	-0.21	-0.56	2.22
tsc_1	1	2.65	0.79	1	2.00	3.00	3.00	4	▁▆▁▇▂	-0.16	-0.40	2.08
tsc_2	1	2.30	0.81	1	2.00	2.00	3.00	4	▂▇▁▆▁	0.19	-0.45	2.11
tsc_3	1	2.67	0.80	1	2.00	3.00	3.00	4	▁▆▁▇▂	-0.20	-0.39	2.10
exp_1	1	5.08	1.42	1	4.00	5.00	6.00	7	▁▁▅▅▇	-0.54	-0.16	2.87
exp_2	1	5.16	1.41	1	4.00	5.00	6.00	7	▁▁▅▃▇	-0.45	-0.44	2.95
exp_3	1	5.16	1.42	1	4.00	5.00	6.00	7	▁▁▅▃▇	-0.53	-0.19	2.92
exp_4	1	5.15	1.46	1	4.00	5.00	6.00	7	▁▂▃▃▇	-0.53	-0.35	2.84
exp_5	1	5.04	1.38	1	4.00	5.00	6.00	7	▁▁▅▅▇	-0.40	-0.28	2.93
exp_6	1	5.01	1.46	1	4.00	5.00	6.00	7	▁▂▅▃▇	-0.42	-0.47	2.75
int_1	1	4.94	1.42	1	4.00	5.00	6.00	7	▁▂▆▅▇	-0.38	-0.34	2.76
int_2	1	4.99	1.39	1	4.00	5.00	6.00	7	▁▂▆▅▇	-0.29	-0.50	2.87
int_3	1	4.73	1.33	1	4.00	5.00	6.00	7	▁▂▇▅▇	-0.24	0.05	2.80
int_4	1	4.89	1.34	1	4.00	5.00	6.00	7	▁▂▇▆▇	-0.29	-0.16	2.90
ben_1	1	4.80	1.36	1	4.00	5.00	6.00	7	▁▂▇▅▇	-0.21	-0.25	2.80
ben_2	1	4.87	1.38	1	4.00	5.00	6.00	7	▁▂▇▅▇	-0.28	-0.25	2.80
ben_3	1	4.89	1.42	1	4.00	5.00	6.00	7	▂▂▆▅▇	-0.35	-0.39	2.74
ben_4	1	4.85	1.32	1	4.00	5.00	6.00	7	▁▂▆▇▇	-0.24	-0.28	2.92
tch_1	1	-176.59	384.37	-999	1.00	3.00	3.00	4	▂▁▁▁▇	-1.67	0.79	2.14
tch_2	1	-197.98	401.39	-999	1.00	3.00	4.00	4	▂▁▁▁▇	-1.49	0.23	2.00
tch_3	1	-264.33	443.31	-999	-999.00	2.00	3.00	4	▃▁▁▁▇	-1.05	-0.89	1.66
tch_4	1	-159.05	368.96	-999	1.00	3.00	3.00	4	▂▁▁▁▇	-1.83	1.37	2.28
tch_5	1	-240.87	430.18	-999	1.00	2.00	3.00	4	▂▁▁▁▇	-1.19	-0.58	1.76
meas_rep	1	1.50	0.50	1	1.00	1.50	2.00	2	▇▁▁▁▇	0.00	-2.00	1.00
age	1	35.54	14.29	16	16.00	35.00	50.00	50	▆▁▅▁▇	-0.34	-1.50	1.37
Experitse	1	5.10	1.27	1	4.17	5.00	6.00	7	▁▁▆▆▇	-0.42	-0.15	3.22
Integrity	1	4.88	1.19	1	4.00	4.75	5.75	7	▁▂▇▇▅	-0.18	-0.02	3.26
Benevolence	1	4.85	1.20	1	4.00	4.75	5.75	7	▁▂▇▇▅	-0.13	-0.28	3.21
TSM	1	2.65	0.57	1	2.25	2.75	3.00	4	▁▂▇▃▂	0.08	0.16	2.92

3.2 METI

3.2.1 Dimensionality (CFA)

First we specified two consecutive threedimensional CFA models

# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                             int_1 + int_2 + int_3 + int_4 +
                             ben_1 + ben_2 + ben_3 + ben_4"

cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 2))

cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4"

cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4"

cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_rbt%>%filter(meas_rep == 1))

## Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
##                 is not positive definite;
##                 use lavInspect(fit, "cov.lv") to investigate.

# define a function which prints the fit
fpf <- function(x){  
  fm_tmp <- fitmeasures(x)
  return(cat(sprintf(
          "χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
           round(fm_tmp[c("chisq")],3), 
                 fm_tmp[c("df")],
           round(fm_tmp[c("cfi")],3),
           round(fm_tmp[c("tli")],3),
           round(fm_tmp[c("rmsea")],3),
           round(fm_tmp[c("srmr")],3),
           round(fm_tmp[c("srmr_between")],3),
           round(fm_tmp[c("srmr_within")],3))))
}

# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)

χ² = 308.08, df = 77, CFI = 0.931, TLI = 0.918, RMSEA = 0.108, SRMR = 0.04, SRMR_between = NA, SRMR_within = NA

fpf(cfa_meti_1)

χ² = 194.143, df = 74, CFI = 0.964, TLI = 0.956, RMSEA = 0.079, SRMR = 0.031, SRMR_between = NA, SRMR_within = NA

cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)

χ² = 297.077, df = 77, CFI = 0.941, TLI = 0.93, RMSEA = 0.105, SRMR = 0.033, SRMR_between = NA, SRMR_within = NA

fpf(cfa_meti_2)

χ² = 196.021, df = 74, CFI = 0.967, TLI = 0.96, RMSEA = 0.08, SRMR = 0.025, SRMR_between = NA, SRMR_within = NA

anova(cfa1d_meti_1, cfa_meti_1)

## Chi-Squared Difference Test
## 
##              Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_1   74 9528.6 9638.6 194.14                                  
## cfa1d_meti_1 77 9636.5 9735.9 308.08     113.94       3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(cfa1d_meti_2, cfa_meti_2)

## Chi-Squared Difference Test
## 
##              Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_2   74 8883.3 8993.3 196.02                                  
## cfa1d_meti_2 77 8978.4 9077.8 297.08     101.06       3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In an next step we ran a two-level CFA …

# onedimensional model
mcfa1d_meti_model <- "level: 1
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4
                    
                    level: 2
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4"



mcfa_meti_model <- "level: 1
                    exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_w =~ int_1 + int_2 + int_3 + int_4 
                    ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
                    
                    level: 2
                    exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_b =~ int_1 + int_2 + int_3 + int_4 
                    ben_b =~ ben_1 + ben_2 + ben_3 + ben_4"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_rbt, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_rbt, cluster = "session")
fpf(mcfa1d_meti)

χ² = 437.579, df = 154, CFI = 0.953, TLI = 0.944, RMSEA = 0.06, SRMR = 0.159, SRMR_between = 0.052, SRMR_within = 0.107

fpf(mcfa_meti)

χ² = 313.519, df = 148, CFI = 0.972, TLI = 0.966, RMSEA = 0.047, SRMR = 0.089, SRMR_between = 0.04, SRMR_within = 0.049

anova(mcfa1d_meti, mcfa_meti)

## Chi-Squared Difference Test
## 
##              Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## mcfa_meti   148 18217 18539 313.52                                  
## mcfa1d_meti 154 18329 18626 437.58     124.06       6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2.2 Reliability (McDonalds \(\omega\))

# First Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("exp")))$est

## [1] 0.9492506

# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("int")))$est

## [1] 0.8795378

# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("ben")))$est

## [1] 0.8819123

# Second Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("exp")))$est

## [1] 0.9504771

# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("int")))$est

## [1] 0.900965

# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("ben")))$est

## [1] 0.9137172

3.3 Topic specific multiplism

3.3.1 Dimensionality (CFA)

First we specified two consecutive onedimensional CFA models

cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
                  tsm_2 ~~ tsm_4"

cfa_tsm_1 <-
  cfa(cfa_tsm_model, 
      data = data_rbt %>% 
        filter(meas_rep == 1),
      missing = "fiml")
fpf(cfa_tsm_1)

χ² = 3.991, df = 4, CFI = 1, TLI = 1, RMSEA = 0, SRMR = 0.03, SRMR_between = NA, SRMR_within = NA

cfa_tsm_2 <-
  cfa(cfa_tsm_model, 
      data = data_rbt %>% 
        filter(meas_rep == 2))
fpf(cfa_tsm_2)

χ² = 5.472, df = 4, CFI = 0.988, TLI = 0.982, RMSEA = 0.038, SRMR = 0.05, SRMR_between = NA, SRMR_within = NA

In an next step, we ran a two-level CFA …

mcfa_tsm_model <- "level: 1
                    tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ tsm_4
                    

                    level: 2
                    tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ 0*tsm_2
tsm_1 ~~ tsm_2"

mcfa_tsm <- cfa(mcfa_tsm_model, data = data_rbt, cluster = "session")
fpf(mcfa_tsm)

χ² = 4.328, df = 3, CFI = 0.995, TLI = 0.98, RMSEA = 0.029, SRMR = 0.092, SRMR_between = 0.071, SRMR_within = 0.02

3.3.2 Reliability (McDonalds \(\omega\))

# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==1) %>% select(starts_with("tsm")))$est

## [1] 0.6962778

# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==2) %>% select(starts_with("tsm")))$est

## [1] 0.6046964

3.4 Treatment check

3.4.1 Dimensionality (CFA)

We specified two consecutive onedimensional CFA models

cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5
                  tch_1 ~~ tch_2"

cfa_tch_model2 <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5"

cfa_tch_1 <- cfa(cfa_tch_model, 
                 data = data_tch_n%>%
                   filter(meas_rep == 1))
fpf(cfa_tch_1)

χ² = 3.982, df = 4, CFI = 1, TLI = 1, RMSEA = 0, SRMR = 0.012, SRMR_between = NA, SRMR_within = NA

cfa_tch_2 <- cfa(cfa_tch_model2, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)

χ² = 8.846, df = 5, CFI = 0.996, TLI = 0.991, RMSEA = 0.066, SRMR = 0.016, SRMR_between = NA, SRMR_within = NA

3.4.2 Reliability (Ordinal McDonalds \(\omega\))

# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 1) %>% select(starts_with("tch")))$est

## [1] 0.8354668

# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est

## [1] 0.9357098

3.5 Table of (M)CFA fit-indices

tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
  column_to_rownames(var = "rownames")%>%
  knitr::kable(., digits = 3)

	1d CFA METI 1	1d CFA METI 2	3d CFA METI 1	3d CFA METI 2	1d MCFA METI	3d MCFA METI	1d CFA TSM 1	1d CFA TSM 2	1d MCFA TSM	1d CFA TCH 1	1d CFA TCH 2
χ²	308.080	297.077	194.143	196.021	437.579	313.519	3.991	5.472	4.328	3.982	8.846
df	77.000	77.000	74.000	74.000	154.000	148.000	4.000	4.000	3.000	4.000	5.000
CFI	0.931	0.941	0.964	0.967	0.953	0.972	1.000	0.988	0.995	1.000	0.996
TLI	0.918	0.930	0.956	0.960	0.944	0.966	1.000	0.982	0.980	1.000	0.991
RMSEA	0.108	0.105	0.079	0.080	0.060	0.047	0.000	0.038	0.029	0.000	0.066
SRMR	0.040	0.033	0.031	0.025	0.159	0.089	0.030	0.050	0.092	0.012	0.016
SRMR_between	NA	NA	NA	NA	0.052	0.040	NA	NA	0.071	NA	NA
SRMR_within	NA	NA	NA	NA	0.107	0.049	NA	NA	0.020	NA	NA
BIC	9735.922	9077.757	9638.632	8993.349	18625.957	18539.351	2421.423	2341.021	4734.532	1709.467	1853.905
AIC	9636.548	8978.383	9528.610	8883.327	18329.002	18216.942	2385.932	2319.726	4645.445	1676.277	1822.088

4 Results of the treatmentcheck

4.1 Plot

res_tch_data <- data_tch%>%
  gather(variable, value, starts_with("tch_"))%>%
  group_by(treat, variable, value)%>%
  mutate(value = as.character(value)) %>% 
  summarize(freq = n())%>%
  ungroup()%>%
  mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
                           treat == "CB" ~ "Colored badges",
                           treat == "CC" ~ "Control condition",
                           T ~ treat),
         value = case_when(value =="-999" ~ "don't know",
                           T ~ value),
         variable = case_when(variable == "tch_1" ~ "item 1", 
                              variable == "tch_2" ~ "item 2", 
                              variable == "tch_3" ~ "item 3", 
                              variable == "tch_4" ~ "item 4", 
                              variable == "tch_5" ~ "item 5"),
         Frequency = freq)

## `summarise()` regrouping output by 'treat', 'variable' (override with `.groups` argument)

res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) + 
  geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
  scale_size_continuous(range = c(3,15)) + 
  scale_color_gradient(low = "grey95", high = "grey65") +
  guides(color=guide_legend(), size = guide_legend()) +
  facet_wrap(~treat) + 
  theme_ipsum_ps() + 
  ggtitle("Results of the treatment check", "Frequency per item and experimental condition") + 
  ylab("") + 
  xlab("")
  
res_tch_plot

res_tch_plot_pub <- ggplot(res_tch_data %>%
                         mutate(treat = case_when(treat == "Greyed out badges" ~ "Grey badges",
                                                  TRUE ~ treat)), aes(variable, value)) + 
  geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
  scale_size_continuous(range = c(3,15)) + 
  scale_color_gradient(low = "grey95", high = "grey65") +
  guides(color=guide_legend(), size = guide_legend()) +
  facet_wrap(~treat) + 
  theme_ipsum_ps() + 
  ylab("") + 
  xlab("")

#ggsave("Fig/res_tch_plot_public.svg", res_tch_plot, width = 11, height = 6)
#ggsave("8_reproducible_documentation_of_analyses/Fig/Fig7.jpg", res_tch_plot_pub, width = 130, height = 50, units = "mm", dpi = 300, scale = 2.4)

4.2 Effect size

res_tch_data_A <- data_tch_n%>%
  filter(treat != "CC")%>%
  filter(tch_1 != -999)

effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)

## 
## Vargha and Delaney A
## 
## A estimate: 0.72609 (medium)

5 Graphical exploration

5.1 Plot Hyp 1

data_scales_lables%>%
  gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
  ggplot(., aes(x = Treatment, y = Value)) + 
  geom_violin(adjust = 1.5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  facet_wrap(~ Variable, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 1)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  ylim(1,7) +
  hrbrthemes::theme_ipsum_ps()

# Export data for combined plot
#data_scales_lables %>% 
#  mutate(Study = "Study 3") %>% 
#  write_csv(., "8_reproducible_documentation_of_analyses/Fig/data_scales_lables_study_3.csv")

5.2 Descriptive Effect Sizes Hyp 1

A_GB_CC <- data_scales_lables%>%
  filter(treat != "CB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_CC_CB <- data_scales_lables%>%
  filter(treat != "GB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_GB_CB <- data_scales_lables%>%
  filter(treat != "CC")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate

	GB	CC
CC	A = 0.56, d = 0.21
CB	A = 0.56, d = 0.2	A = 0.49, d = -0.02

5.3 Hyp 2/3

plot_hyp2_1 <- ggplot(data_scales_lables, 
                      aes(x=`Topic specific multiplism`, y = Integrity)) + 
  geom_jitter() +
  facet_wrap(~ Treatment, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 2/3)",
       subtitle = "Jitter plot per treatment") +
  hrbrthemes::theme_ipsum_ps()


plot_hyp2_1 + stat_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plot_hyp2_1 + stat_smooth(method = "lm")

## `geom_smooth()` using formula 'y ~ x'

5.4 Descriptive Effect Sizes Hyp 3/4

Spearman and Kendall correlations:

r_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

	GB	CC	CB
\(r(integrity, multiplism)\)	-0.05	-0.13	0.05
\(\tau(integrity, multiplism)\)	-0.04	-0.07	0

5.5 Hyp 4

data_scales_lables%>%
  mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
                               treat == "CB" ~ "Colored badges",
                               treat == "CC" ~ "Control condition",
                               T ~ "treat"))%>%
  ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
  geom_violin(adjust = 1.5, alpha = .5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  labs(title = "Graphical overview (Hyp 4)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  xlab("") +
  ylim(1,4) +
  hrbrthemes::theme_ipsum_ps()

fig4 <- data_scales_lables%>%
 mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
                              treat == "CB" ~ "Colored badges",
                              treat == "CC" ~ "Control condition",
                              T ~ "treat"))%>%
 ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
 geom_violin(adjust = 1.5, alpha = .5) +
 stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
 coord_flip() +
 xlab("") +
 ylim(1,4) +
 hrbrthemes::theme_ipsum_ps()

5.6 Descriptive Effect Sizes Hyp 4

A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)


A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "GB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)


A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CC")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)

	GB	CC
CC	A = 0.49, d = -0.02
CB	A = 0.46, d = -0.15	A = 0.46, d = -0.13

6 Inference statistics (Hyp 1)

6.1 Description of the variables

All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:

GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open Practices
CC means, that the participants of the study could not learn if or if not the authors of the study used Open Practices
CBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open Practices

Finally, the variable session identified the study participants.

If we look descriptively at these variables:

data_scales_wide%>%
  my_skim(.)

Data summary
Name	Piped data
Number of rows	257
Number of columns	13
_______________________
Column type frequency:
character	1
numeric	12
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
session	0	1	64	64	0	257	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist	skewn	kurto	maxabszvalue
Benevolence_CB	96	0.63	4.94	1.17	1.75	4.00	5.00	5.75	7	▁▃▇▆▆	0.03	-0.66	2.74
Benevolence_CC	87	0.66	4.96	1.22	1.25	4.00	5.00	6.00	7	▁▂▇▇▇	-0.34	-0.04	3.05
Benevolence_GB	74	0.71	4.68	1.20	1.00	4.00	4.75	5.50	7	▁▂▇▇▅	-0.06	-0.21	3.07
Experitse_CB	96	0.63	5.11	1.28	1.00	4.00	5.00	6.00	7	▁▂▆▆▇	-0.33	-0.38	3.21
Experitse_CC	87	0.66	5.21	1.23	1.00	4.33	5.33	6.00	7	▁▁▅▆▇	-0.63	0.40	3.43
Experitse_GB	74	0.71	4.99	1.31	1.00	4.00	5.00	6.00	7	▁▂▇▇▇	-0.31	-0.37	3.04
Integrity_CB	96	0.63	4.96	1.15	1.25	4.00	5.00	5.75	7	▁▂▆▇▅	-0.11	-0.21	3.24
Integrity_CC	87	0.66	4.97	1.19	1.00	4.00	5.00	6.00	7	▁▁▇▇▆	-0.30	0.21	3.33
Integrity_GB	74	0.71	4.73	1.22	1.00	4.00	4.50	5.75	7	▁▂▇▆▅	-0.11	-0.13	3.06
TSM_CB	96	0.63	2.60	0.58	1.00	2.25	2.50	3.00	4	▁▃▇▃▂	0.09	0.06	2.74
TSM_CC	87	0.66	2.67	0.57	1.00	2.25	2.75	3.00	4	▁▂▇▃▂	0.04	0.20	2.93
TSM_GB	74	0.71	2.68	0.55	1.00	2.25	2.75	3.00	4	▁▂▇▅▂	0.14	0.13	3.08

6.2 Investigating the missingness

6.2.1 Missingness per Variable

library(naniar)

## 
## Attaching package: 'naniar'

## The following object is masked from 'package:skimr':
## 
##     n_complete

visdat::vis_miss(data_scales_wide%>%select(-session)) + 
  ggtitle("Missingness per Variable") + 
  theme(plot.margin = margin(1, 2, 1, 1, "cm"))

6.2.2 Marginal distributions `Integrity_GB`

library(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.3 Marginal distributions `Integrity_CC`

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.4 Marginal distributions `Integrity_CB`

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.5 Marginal distributions `TSM_GB`

ggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.6 Marginal distributions `TSM_CC`

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.7 Marginal distributions `TSM_CB`

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.8 Cohen’s d of missing/not-missing per variable

d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
  round(.,4)

knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")

Boxplot of Cohen’s d of missing/not-missing per variable
	Integrity_CB	Integrity_CC	Integrity_GB	TSM_CB	TSM_CC	TSM_GB
Integrity_CB	NaN	-0.1487	0.1183	NaN	-0.0571	0.0617
Integrity_CC	-0.0627	NaN	-0.1183	-0.1243	NaN	-0.0617
Integrity_GB	0.0627	0.1487	NaN	0.1243	0.0571	NaN
TSM_CB	NaN	-0.1487	0.1183	NaN	-0.0571	0.0617
TSM_CC	-0.0627	NaN	-0.1183	-0.1243	NaN	-0.0617
TSM_GB	0.0627	0.1487	NaN	0.1243	0.0571	NaN

boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")

6.3 Imputation

M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
            m = M,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

6.4 Check of first 10 imputation

out_first10 <-   mice(data = data_scales_wide%>%select(-session),
            m = 10,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

densityplot(out_first10)

plot(out_first10)

6.5 Parameter and BF estimation

# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3) 
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3) 


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance


# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))* 
  sum(diag(covbetween_hyp1 %*% 
             MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)


# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp1 <- bain(estimates_hyp1,
               "Integrity_GB<Integrity_CC<Integrity_CB;
                Integrity_GB=Integrity_CC=Integrity_CB;
                Integrity_GB<Integrity_CC=Integrity_CB",
                n = neff_hyp1, Sigma=covariance_hyp1,    
                group_parameters=3, joint_parameters = 0)
print(results_hyp1)

## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.516 0.161 3.199  5.543  0.135 0.130
## H2 0.476 0.237 2.012  2.012  0.085 0.082
## H3 4.011 0.217 18.473 18.473 0.780 0.748
## Hu                                 0.041
## 
## Hypotheses:
##   H1: Integrity_GB<Integrity_CC<Integrity_CB
##   H2: Integrity_GB=Integrity_CC=Integrity_CB
##   H3: Integrity_GB<Integrity_CC=Integrity_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.

summary(results_hyp1)

##      Parameter        n Estimate       lb       ub
## 1 Integrity_GB 172.6037 4.746086 4.581740 4.910432
## 2 Integrity_CC 172.6037 4.958640 4.794640 5.122640
## 3 Integrity_CB 172.6037 4.965926 4.802789 5.129063

results_hyp1$BFmatrix

##           H1       H2        H3
## H1 1.0000000 1.589591 0.1731532
## H2 0.6290927 1.000000 0.1089294
## H3 5.7752334 9.180258 1.0000000

7 Inference statistics (Hyp 2)

7.1 Parameter estimation

path_mod <- "Integrity_GB ~ TSM_GB
             Integrity_CC ~ TSM_CC             
             Integrity_CB ~ TSM_CB"

# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 path_fit <- sem(path_mod, 
                 data = mice::complete(out, i), 
                 std.lv = T) # estimate the path coefficients 
 best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
                     filter(op == "~")%>%
                     pull(std.all) 
 covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
                  1/M * lavInspect(path_fit, 
                                   "vcov.std.all")[c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB"),
                                                   c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB")]
}

# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)

## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB 
##         -0.03         -0.13         -0.09

7.2 Visual path model

\usetikzlibrary{arrows,positioning}
\usetikzlibrary{calc}
\begin{tikzpicture}[auto,>=latex,align=center,
latent/.style={circle,draw, thick,inner sep=0pt,minimum size=10mm},
manifest/.style={rectangle,draw, thick,inner sep=2pt,minimum size=15mm},
manifestcov/.style={rectangle,draw, thick,inner sep=1pt,minimum size=8mm},
mean/.style={regular polygon,regular polygon sides=3,draw,thick,inner sep=0pt,minimum size=10mm},
paths/.style={->, thick, >=latex, font = \footnotesize\sffamily, fill = white},
variance/.style={<->, thin, >=latex, bend left=270, looseness=2.5, font = \footnotesize},
covariance/.style={<->, thin, >=latex, bend left=290, looseness=.5, font = \footnotesize},
dotsdots/.style={circle, draw, fill = black, minimum size = 1mm, maximum size = 1mm}
]

%% AVs
\node [manifest] (AV1) at (0, -2)       {$int_{gb}$};
\node [manifest] (AV2) [right = of AV1] {$int_{cc}$};
\node [manifest] (AV3) [right = of AV2] {$int_{cb}$};


%% UVs
\node [manifest] (UV1) at (0, 2)        {$tsm_{gb}$};    
\node [manifest] (UV2) [right =of UV1]  {$tsm_{cc}$};    
\node [manifest] (UV3) [right =of UV2]  {$tsm_{cb}$};    

%% Paths
\draw [paths] (UV1) to node [midway]{-.03} (AV1);
\draw [paths] (UV2) to node [midway]{-.13} (AV2);
\draw [paths] (UV3) to node [midway]{-.09} (AV3);

%\draw [covariance] (UV2) to node {} (UV1);
\end
{tikzpicture}

covbetween_hyp2 <- cov(best_hyp2) # between covariance matrix
covariance_hyp2 <- covwithin_hyp2 + (1+1/M)*covbetween_hyp2 # total variance

# Determine the effective and real sample sizes ######
samp_hyp2 <- nrow(data_scales_wide) # real sample size
nucom_hyp2 <- samp_hyp2-length(estimates_hyp2)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp2 <- (1+1/M)*(1/length(estimates_hyp2))* 
  sum(diag(covbetween_hyp2 %*% 
             MASS::ginv(covariance_hyp2))) # ... (43)
nuold_hyp2 <- (M-1)/(lam_hyp2^2) # ... (44)
nuobs_hyp2 <- (nucom_hyp2+1)/(nucom_hyp2+3)*nucom_hyp2*(1-lam_hyp2) # ... (46)
nu_hyp2 <- nuold_hyp2*nuobs_hyp2/(nuold_hyp2+nuobs_hyp2) # ... (47)
fracmis_hyp2 <- (nu_hyp2+1)/(nu_hyp2+3)*lam_hyp2 + 2/(nu_hyp2+3) # ... (48)
neff_hyp2 <- samp_hyp2-samp_hyp2*fracmis_hyp2 # = 114 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp2 <- list(covariance_hyp2)

7.3 Bayes Factor estimation (Hyp 2)

set.seed(6346)
results_hyp2 <- bain(estimates_hyp2, 
                     "Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
                     Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = F)
print(results_hyp2)

## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.570  0.163 3.495  6.804  0.094 0.092
## H2 22.961 0.685 33.497 33.497 0.906 0.882
## Hu                                  0.026
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
##   H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.

results_hyp2$BFmatrix

##          H1        H2
## H1 1.000000 0.1043274
## H2 9.585207 1.0000000

8 Inference statistics (Hyp 3)

set.seed(6346)
results_hyp3 <- bain(estimates_hyp2, 
                    "Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
                     (Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
                     Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = T)

print(results_hyp3)

## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.050  0.169 0.295  0.258  0.009 0.009
## H2 0.259  0.346 0.750  0.663  0.023 0.023
## H3 10.821 0.347 31.201 31.201 0.968 0.938
## Hu                                  0.030
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
##   H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
##   H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.

results_hyp3$BFmatrix

##           H1         H2          H3
## H1   1.00000  0.3938263 0.009469693
## H2   2.53919  1.0000000 0.024045352
## H3 105.60005 41.5880783 1.000000000

9 Inference statistics (Hyp 4)

# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance


# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))* 
  sum(diag(covbetween_hyp4 %*% 
             MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)


# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp4 <- bain(estimates_hyp4,
               "TSM_GB>TSM_CC>TSM_CB;
                TSM_GB=TSM_CC=TSM_CB;
                (TSM_GB,TSM_CC)>TSM_CB",
                n = neff_hyp4, Sigma=covariance_hyp4,    
                group_parameters=3, joint_parameters = 0)

print(results_hyp4)

## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.619  0.176 3.527  7.642  0.125 0.120
## H2 17.926 0.799 22.432 22.432 0.792 0.765
## H3 0.826  0.350 2.359  8.796  0.083 0.080
## Hu                                  0.034
## 
## Hypotheses:
##   H1: TSM_GB>TSM_CC>TSM_CB
##   H2: TSM_GB=TSM_CC=TSM_CB
##   H3: (TSM_GB,TSM_CC)>TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.

results_hyp4$BFmatrix

##           H1        H2       H3
## H1 1.0000000 0.1572512 1.495238
## H2 6.3592510 1.0000000 9.508594
## H3 0.6687899 0.1051680 1.000000

Reproducible Documentation of Analyses

Public sample